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Abstract. In this paper we conduct a systematic study of the granularity of the 
initial state of hot and dense QCD matter produced in ultra-relativistic heavy-ion 
collisions and its influence on bulk observables like particle yields, itit spectra and 
elliptic flow. For our investigation we use a hybrid transport model, based on 
(3-|-l)d hydrodynamics and a microscopic Boltzmann transport approach. The initial 
conditions are generated by a non-equilibrium hadronic transport approach and the 
size of their fluctuations can be adjusted by defining a Gaussian smoothing parameter 
a. The dependence of the hydrodynamic evolution on the choices of a and tstart is 
explored by means of a Gaussian emulator. To generate particle yields and elliptic flow 
that are compatible with experimental data the initial state parameters are constrained 
to be (7 = 1 fm and tstart = 0.5 fm. In addition, the influence of changes in the equation 
of state is studied and the results of our event-by-event calculations are compared to a 
calculation with averaged initial conditions. We conclude that even though the initial 
state parameters can be constrained by yields and elliptic fiow, the granularity needs 
to be constrained by other correlation and fluctuation observables. 
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1. Introduction 

The study of the properties of strongly interacting matter at very high temperatures 
and/or densities is the goal of the experimental heavy ion program at the Relativistic 
Heavy Ion Collider (RHIC) at Brookhaven National Laboratory (BNL) and at the Large 
Hadron Collider (LHC) at the European Organization for Nuclear Research (CERN) 
in Geneva. In order to extract the relevant QCD properties, such as the viscosity and 
other transport coefficients, from the final state particle distributions measured by the 
experiments, effective dynamical approaches are needed that connect theoretical input 
with experimental data. 

So called hybrid approaches that apply ideal relativistic fiuid dynamics for the 
hot and dense stage of the collision, when the produced matter is close to thermal 
equilibrium, and use microscopic transport approaches for the initial and/or final non- 
equilibrium stages have proven to be the most successful approaches for the description 
of relativistic heavy-ion collisions [H El [3l HJ El E]. Unfortunately, as the models gain 
complexity and sophistication, the number of parameters used to encode the relevant 
physics increases. All of these parameters need to be constrained or determined before 
conclusions about the properties of hot and dense QCD matter can be drawn - this 
constitutes a non-trivial multi-parameter optimization problem. 

The Models and Data Analysis Initiative (MADAI) collaboration has been formed 
with the aim of developing novel techniques based on statistical science to address 
multiparameter estimation problems one encounters in the application of complex 
models to the extraction of knowledge from large multi-dimensional experimental data 
sets. Here, we will focus on ultra-relativistic heavy-ion collisions, in particular on the 
problem of constraining the initial conditions of such collisions via a comparison of bulk 
observables calculated in our model approach to data. 

Recently, there has been a rising interest in quantifying the initial state fiuctuations 
in ultra-relativistic heavy-ion collisions due to their importance for multi-particle 
correlation measurements in particular regarding elliptic fiow, triangular fiow and the 
so-called "ridge" correlations [3 [HI El [TOl [IH [I2] . 

In this paper a systematic study of the initial state granularity and its infiuence 
on bulk observables such as particle yields, spectra and elliptic fiow is presented. 
In Section [2] the employed hybrid transport approach is introduced and in Section El 
the statistical methods are discussed. In Section H] we describe how particle yields 
and elliptic fiow can be used to constrain the initial state parameters. The infiuence 
of changes in the equation of state and a comparison to averaged initial conditions is 
contained in Section El In the last Section El the conclusions from this systematic study 
are summarized. 
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2. The Hybrid Approach 

The present hybrid model used to simulate the dynamics of Au+Au collisions at 
^/suN = 200 GeV is based on UrQMD [131 El with an intermediate ideal hydrodynamic 
evolution for the hot and dense stage of the reaction [5J. UrQMD is a string/hadronic 
transport model which simulates multiple interactions of ingoing and newly produced 
particles, the excitation and fragmentation of color strings flEi [T6] and the formation 
and decay of hadronic resonances. To mimic experimental conditions as realistically as 
possible the non-equilibrium dynamics in the initial and the final state are taken into 
account on an event-by-event basis (see e.g. discussions in [T7 1 [T8l \T9\ [20 1 [211 122] ). 

In UrQMD, the incoming nuclei are initialized according to Woods-Saxon profiles 
and the initial nucleon-nucleon scatterings and non-equilibrium dynamics proceed 
according to the Boltzmann equation. After the two nuclei have passed through 
each other local thermal equilibrium is assumed to make the transition to the ideal 
hydrodynamic description. The time chosen for the mapping of the particle degrees of 
freedom to the thermodynamic fields is called tgtart- This is one of the parameters that 
we need to adjust in this calculation (default is tstart = 0.5 fm). The point particles are 
represented by three-dimensional Gaussian distributions in the following way 

3 

e(x,i/,z)=^-J -,Epe.p — (1) 

to obtain energy, momentum and net baryon density distributions that are smooth 
enough for the hydrodynamic evolution. Here, e is the energy density at position (x, y, z) 
that a particle with energy Ep at position {xp,yp,Zp) contributes. The Gaussians are 
Lorentz contracted in z-direction by 7^ to account for the large longitudinal velocities. 
Only the matter at midrapidity \y\ < 2 is assumed to be locally equilibrated, whereas the 
other hadrons are treated in the hadronic cascade. The second parameter that influences 
the smearing and therefore the granularity of the initial condition is the width of the 
Gaussian a. This parameter will be varied between a = 0.8 fm, the lower limit to keep 
the hydrodynamic code numerically stable and a = 2 fm which leads to very smooth 
profiles as shown in ^23j . The default choice is cr = 1 fm which corresponds to the typical 
size of a nucleon. 

In this manner, the initial conditions for the hydrodynamic calculation include 
fluctuations from the early non-equilibrium dynamics such as finite velocity profiles and 
peaks in the energy density because of fluctuations in the energy deposition. Starting 
from these (single event) initial conditions a full (3-1-1) dimensional ideal hydrodynamic 
evolution is performed using the SHASTA algorithm [241 EH]. The hydrodynamic 
evolution is stopped if a certain transition criterion is fulfilled. In [26] we have explored 
a freeze-out procedure to account for the large time dilatation that occurs for fluid 
elements at large rapidities. To mimic an iso-r hypersurface we freeze out full transverse 
slices, of thickness Az = 0.2fm, whenever all cells of that slice fulfill the freeze-out 
criterion (e ~ 713 MeV/fm^). On these slices particles are produced according to the 
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Cooper-Frye formula and the hadronic rescattering and resonance decays are taken into 
account in the hadronic cascade (UrQMD). 

The event-by-event calculation provides the full final phase-space distribution of 
the hadrons that are also measured in experiments and therefore allows for detailed 
comparisons of many observables at the same time. In the present study we concentrate 
on the question how the initial state granularity can be constrained by bulk observables. 
In the following section, we shall demonstrate how a more sophisticated statistical 
analysis helps to get a good handle on this multi-parameter fit problem. 

Other important parameters that influence the results are finite shear and bulk 
viscosity, the freeze-out criterion and the equation of state that is employed for the 
calculation. Some of these parameters are not implemented in our current model (e.g. 
the shear and bulk viscosities) and therefore an extensive study of these additional 
parameter dependencies is left to a future publication. 

3. Gaussian Process Emulators for Computer Models 

A thorough exploration of the dependence of UrQMD+hydrodynamics predictions of 
transverse mass spectra for midrapidity pions and kaons on such model parameters 
as a and t^^_^^^ could entail thousands of code evaluations taking each ~ 3 CPU- 
hours. As an alternative, we construct a model emulator [271 [28] — a statistical 
model of UrQMD+hydrodynamics which, after "training" with the model results at 
a designed finite set of locations in its parameter space, can predict rapidly what 
the hybrid approach predictions would be at other untried parameter vectors, with 
an attendant measure of uncertainty. Unlike simpler interpolation schemes, which 
typically provide only a single estimated value at each point in the space, the emulator 
generates an entire joint probability distribution representing what is known about 
the unobserved computer model outputs, supporting the generation of predictions of 
arbitrary functions of the computer model outputs with measures of their uncertainty. 
Numerical implementations of the emulator models are so fast that it is entirely practical 
to simulate thousands or even millions of predicted model outputs in minutes of 
computer time. This makes possible a broad range of interesting analysis of the output 
of computer codes which would otherwise require an unacceptably large computational 
effort (see [291 EQl |3T] , for example). 

For our emulator the computer model's output at all possible input vectors {9} is 
modeled as a Gaussian stochastic process or random field, indexed by 6, expressing initial 
uncertainty about possible model results; several random draws from such a distribution 
are illustrated for a one-dimensional parameter space in Fig. [Deleft). Emulator training 
is accomplished by evaluating the conditional probability distribution of model outputs 
at all input vectors {9}, given observations of the actual computer model output at a 
finite collection V of selected design input vectors {9i G V}. The emulator reproduces 
the model output perfectly at the design points, and offers predictive distributions for 
model output at other points which are more (or less) variable for points that are closer to 
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(or further from) the design points, respectively. Although other distributions could be 
used for emulators, Gaussian processes (or GPs) are particularly convenient because the 
required conditional distributions can be computed very efficiently with simple matrix 
algebra. 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 



Figure 1. Left: unconditioned draws ft-om a Gaussian proeess. Riglit: draws ft-om 
the same process after conditioning on 7 training points (black circles) from a simple 
model. The grey band is a 95% confidence interval. Note how the uncertainty grows 
away from the training points. 

The detailed behaviour of the emulator is determined by the mean and covariance 
structure chosen for the GP and by the number and location of the training points. 
After removing simple linear effects with multiple regression and a possible nonlinear 
transformation of the computer model's inputs and outputs, if necessary, the computer 
model may be assumed to be locally smooth (so small changes in 6 lead to small 
changes in the output) and homogeneous (so the magnitude of predictive uncertainty 
is relatively constant). Under these assumptions the conventional mean- zero isotropic 
Gaussian random fields used in geostatistics [321 [33] work well; we used the commonly- 
recommended power exponential family [Ml |35] , with power close to its upper limit of 
two to ensure smoothness, and with correlation lengths fit to the data using maximum 
likelihood methods (emulator dependence on this parameter is illustrated in Fig. |2]). A 
modest number of training points, perhaps 10-15 points per dimension in the parameter 
space, evenly distributed over the parameter range of interest will suffice. We used a 
Latin hypercube design, a space-filling approach which has proven very successful for 
all-purpose designs of computer experiment runs because it can "fill" the design space 
with very few points (see [35] or [361 Chap. 5]). 

To emulate functional output, when the computer model output for each input 
vector 6 is not a single numerical quantity but a function of one or more variables. 
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Figure 2. Varying the correlation length 9 changes the emulator dramatically. 
Solid blue lines represent emulator means, while pointwise 95% (pointwise) predictive 
intervals are shown in grey. Left shows over-fitting with 6 too small; right shows 
over-smoothing with 9 too large; middle shows optimal value oi 9 = 0.68, chosen by 
likelihood maximization. 



such as rapidity distributions or particle spectra which vary as a function of centrahty, 
beam energy, etc., we expand those functional outputs in an orthogonal basis (the 
most efficient choice is to use Principal Components), then construct independent 
Gaussian emulators as above for each component separately. This is only marginally 
more computationally demanding than the emulation of univariate computer model 
output. 

In summary: emulators are computationally-trivial statistical approximations to 
functions that open exciting new doors for the statistical exploration of computationally- 
expensive models. 

4. Initial State Variation 

Let us now apply the emulator based on a Gaussian regression process to determine the 
initial state parameters from particle yields and spectra. By varying the Gaussian width 
that has been introduced in Section |2] the initial state granularity can be changed while 
keeping everything else constant. Fig. [3] (left) shows the result for transverse momentum 
spectra for pions and kaons in central Au-|-Au collisions at the highest RHIC energy. 
The slope of the spectra stays constant while the overall height which corresponds to 
the respective particle yield increases with increasing a. This can be traced back to an 
increase in the total entropy of the initial state that is larger, if the particles are smeared 
out over a larger phase space volume. 

From earlier studies in [5] at lower energies it is known that the yields are also 
affected by the choice of the starting time. There is the possibility to compensate a 
change in a by adjusting tstart at the same time. Please note, that the starting time 
is here not only a way to normalize the initial state distributions as in hydrodynamic 
calculations with smooth initial conditions, but it is also the transition time from the 
hadronic transport approach to ideal hydrodynamics. By using the emulator to explore 
the two-dimensional parameter space in a and tstart, we can find parameter combinations 
that lead to the same pion yield/entropy in the system. The following parameter pairs 
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Figure 3. Transverse mass spectra for tt~ and at midrapidity {\y\ < 0.5) in central 
(6 < 3.4 fm) Au+Au collisions at .^/snn = 200 GeV from the hybrid approach for 
different values of a with a fixed starting time (left) and a varied starting time (right) 
represented by lines in comparison to experimental data that is shown as symbols 



can be identified and are listed in table [TJ 



a [fm] 


0.8 


1.0 


1.5 


2.0 


tstart [fm] 


0.3 


0.5 


1.2 


1.8 



Table 1. Combinations of a and tstart that lead to the same pion yield in central 
{b < 3.4 fm) Au+Au collisions at ^snn = 200 GeV. These values have been used in 
all the calculations that are dubbed with 'tstart varied'. 



Fig. m shows the emulated mean number of pions produced in one collision as a 
function of the two initial state parameters tstart and a. A strong linear correlation 
between the parameters can be observed. The emulator was trained on 30 observations 
of the model distributed somewhat uniformly in the 2d parameter space. The variance of 
the emulator (Fig. E]) in the region where the simulation data was collected a e (0.2, 3.0) 
fm and tstart ^ (0.8, 3) fm is small, so we can be confident in the shape of the surface. 
The variance increases for a < 0.8 fm which is a region of pure extrapolation, this is the 
expected behavior. 

The linear correlation between start time and kernel width implies that more 
smearing in the initial state can be compensated for by prolonged evolution before 
the switch to hydrodynamical evolution. 

The result for the pion and kaon transverse momentum spectra for the four 
parameter combinations are shown in Fig. [3] (right). The particle yields are now kept 
constant, but the slope of the spectra is sensitive to the starting time. The earlier the 
starting time, the longer the duration of the hydrodynamic evolution, therefore there is 
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more time to develop larger radial flow which leads to flatter transverse mass spectra. 
The difference in slope becomes clearly visible at higher my — rriQ values larger than 1 
GeV. 



o 




tstart (fm/c) 



Figure 4. Emulated number of pions at midrapidity {\y\ < 0.5) for central (6 < 3.4 
fm) Au+Au collisions at ^snn — 200 GeV in the two-dimensional parameter space of 
a and tstart- 

Let us now explore the influence of the different initial state parameters on the 
integrated elliptic flow in mid-central (6 = 5 — 9 fm) Au+Au collisions at a/snn = 200 
GeV. In Fig. |6]the averaged value of the second coefficient of the Fourier decomposition 
of the azimuthal distribution of charged particles in momentum space is shown as a 
function of the granularity a. The filled circles represent the results for fixed starting 
time tstart = 0.5 fm and the open squares show the result for the identified combinations 
in Table [TJ Without changing the starting time the results for different granularities 
are compatible to the experimental data, that is indicated by the lines, even though the 
particle numbers are very different. The starting time can be constrained by looking at 
the spectra and elliptic flow to be less than 1 fm since later starting times do not allow 
for enough flow development during the hydrodynamic evolution. Just by looking at 
basic bulk observables like yields, spectra and integrated elliptic flow, we can constrain 
the initial state parameters within this approach to be around a = 1 fm and tstart = 0.5 
fm which we will use as default in the following Section. 
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Figure 5. Emulated number of pions at midrapidity {\y\ < 0.5) for central (6 < 3.4 
fm) Au+Au collisions at ^snn — 200 GeV. Left: Transect at fixed tstart and right: at 
fixed a. 
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Figure 6. The averaged value of elliptic flow of charged particles at midrapidity 
(|y| < 0.5) for mid-central (6 = 5 — 9 fm) Au+Au collisions at ^snn — 200 GeV 
as a function of a with fixed (full circles) and varied istart (open squares) compared 
to experimental data represented by black lines (the grey shaded regions indicate the 
error bars) gOl gH |42] . 



5. Fluctuating vs. Averaged Initial Conditions 



Even though, we have shown in the last Section that the initial state parameters can be 
constrained by bulk observables, this does not necessarily constrain the granularity of the 
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initial state. Previous studies at lower energies within the same hybrid approach have 
shown that there is some insensitivity to the event-by-event fluctuations of the initial 
state jl3]. In contrast to the full event-by-event setup we have discussed so far, one can 
also look at calculations from averaged initial conditions. The default parameters are 
chosen to generate the initial state and one averages over 100 initial states from UrQMD 
to feed a smooth profile in the hydrodynamic calculation. This provides a different way 
to tune the initial state granularity between fluctuating (1 UrQMD event) to averaged 
smooth initial conditions (100 UrQMD events). 




Figure 7. Transverse mass spectra for 7r~ and (left) and elliptic flow as a function 
of transverse momentum for pions (right) at midrapidity {\y\ < 0.5) in central/mid- 
central (5 < 3.4 fm/6 = 5 — 9 fm) Au+Au collisions at ^snn ~ 200 GeV from the 
hybrid approach for fluctuating initial conditions with two different EoS and averaged 
initial conditions represented by lines in comparison to experimental data that is shown 
as symbols [33 [Ml ESI |44] . 

In Fig. [7] results for transverse mass spectra (left) and elliptic flow of pions as a 
function of transverse momentum (right) are shown for the fluctuating event-by-event 
setup and the averaged initial conditions. The results in this case are not affected since 
the V2 analysis has been performed with respect to the reaction plane that is given by the 
coordinate system. Using the standard event plane method the elliptic flow results would 
be a little higher in the fluctuating has been shown in |15]. Furthermore, a softer 

equation of state based on a chiral hadronic Lagrangian that is coupled to the Polyakov 
loop including chiral symmetry restoration and a deconfinement phase transition that 
reproduces ground state properties and results from lattice QCD (DE-EoS) [iGj [ITJ \M\ 
has been applied to give an example of other 'parameters' that may affect the results. 
As shown in Fig. [7] (left) the slope of the transverse mass spectra is smaller using the 
DE-EoS and leads to a better agreement with experimental data, while the elliptic flow 
result (right) is very similar. For the calculations with the other equation of state the 
freeze-out transition needs to be adjusted as well and a higher energy density criterion 
has been applied. A more detailed study of the freeze-out procedure is left for a future 
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Figure 8. The yields (left) and the mean transverse mass (right) at midrapidity 
{\y\ < 0.5) for four different particle species (tt", K~^, P, A) are calculated in central 
{b < 3.4 fm) Au+Au collisions at ^snn — 200 GeV. The hybrid approach for 
different values of a with a fixed starting time (full symbols) and a varied starting 
time (dashed lines) are compared to experimental data indicated by full (STAR) and 
dotted (BRAHMS) lines [49l [50l |3Zl [51] . In addition, the crosses represent the event-by- 
event calculation with the equation of state including a deconfinement phase transition 
(DE-EoS) and results for averaged initial conditions are shown as stars. 

To wrap up the results of the present systematic study on initial state granularity 
and bulk observables in heavy ion reactions, in Fig. [8] results for yields and mean 
transverse momentum of pions, protons, kaons and A' s are shown for all the 10 
different cases that we have discussed in this paper so far. The full symbols represent 
the results where only a has been changed while tstart is fixed to be 0.5 fm. The 
multiplicities are proportional to the smearing parameter, whereas the radial flow is 
not as sensitive (with the exception of the A's). By varying the starting time according 
to the result of the emulator (lines), we confirm that the yields can be kept constant but 
the mean transverse mass decreases significantly because of the later transition to the 
hydrodynamic evolution. One needs to start early enough to allow to develop radial flow 
and elliptic flow. From these two calculations the default parameters for a and tstart are 
constrained. Still one might vary the granularity by averaging over initial states instead 
of final states and the corresponding results are shown as stars. Bulk observables do 
not provide a good handle on the granularity since the results are very similar in the 
fluctuating and the averaged case. One needs to calculate fluctuation or correlation 
observables to really constrain the amount of initial state fluctuations. The equation 
of state including a deconfinement phase transition reduces the mean transverse masses 
to achieve a better agreement with the experimental data. We are not able to present 
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a complete set of best fit parameters yet, since the freeze-out procedure needs to be 
improved before the statistical analysis is extended to include more parameters and 
more experimental data sets. 

6. Summary and Conclusion 

We have explored the capabilities of a novel statistical analysis utilizing Gaussian 
process emulators to address ambiguities in the initial state parameters used in a hybrid 
approach for the dynamical evolution of relativistic heavy ion reactions. The smearing 
parameter a and the starting time of the hydrodynamic evolution are found to be 
correlated. By imposing a constant entropy constraint one finds viable combinations of 
these two parameters, which agree well with experimental data. 

The comparison of the fluctuating event-by-event setup to averaged initial 
conditions using the same initial state parameters shows that bulk observables are 
insensitive to the initial state granularity. This insensitivity can be exploited in order to 
constrain quantities like e.g. the shear viscosity of QCD matter from averaged elliptic 
flow results. We note that significant sensitivities to the equation of state remain which 
can be utilized to constrain this input parameter. Our results show that state of the art 
statistical techniques for executing multi-parameter fits are required to reliably quantify 
properties of interest of hot and dense QCD matter. These techniques are now being 
developed and will lead us to achieve a better understanding of sensitivities on different 
parameter sets. 

Acknowledgements 

We are grateful to the Open Science Grid for the computing resources. The authors 
thank Dirk Rischke for providing the 1 fluid hydrodynamics code. H.P. acknowledges 
a Feodor Lynen fellowship of the Alexander von Humboldt foundation. This work 
was supported in part by U.S. department of Energy grant DE-FG02-05ER41367, NSF 
grants PHY-09-41373 and DMS-07-57549 and NASA grant NNX09AK60G. 

References 

[1] S. A. Bass and A. Dumitru, Phys. Rev. C 61, 064909 (2000) | arXiv:nucl-th/0001033] . 
[2] D. Teaney, J. Lauret and E. V. Shuryak, arXiv:nucl-th/0110037l 

[3] T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey and Y. Nara, Phys. Lett. B 636, 299 (2006) 

;arXiv:nucl-th/0511046. 
[4] C. Nonaka and S. A. Bass, Phys. Rev. C 75, 014902 (2007) [arXiv:nucl-th/0607018] . 
[5] H. Petersen, J. Steinheimer, G. Burau, M. Bleicher and H. Stocker, Phys. Rev. C 78, 044901 

(2008) [arXiv:08n6T 695 [nucl-th]]. 
[6] K. Werner, I. Karpenko, T. Pierog, M. Bleicher and K. Mikhailov, Phys. Rev. C 82, 044904 (2010) 

[arXiv: 1004.0805 [nucl-th]]. 
[7] H. Agakishiev et al, a rXiv: 1010.06901 [nucl-ex]. 

[8] B. H. Alver, C. Gombeaud, M. Luzum and J. Y. Olhtrault, Phys. Rev. C 82, 034913 (2010) 
,arXiv: 1007.546 9 [nucl-th]]. 



Initial State Granularity 



13 



[9] B. Schenke, S. Jeon and C. Gale, arXiv: 1009:3241 [hep-ph]. 
[10] H. Petersen, G. Y. Qin, S. A. Bass and B. Muller, Phys. Rev. C 82, 041901 (2010) |arXiv: 1008.06251 
[nucl-th]]. 

[11] G. Y. Qin, H. Petersen, S. A. Bass and B. Muller. [arXiv:1009.1847 [nucl-th]. 
[12] A. Mocsy and P. Sorensen. larXiv:10 08.3381 [hep-ph]. 

[13] S. A. Bass et ai, Prog. Part. Nucl. Phys. 41, 255 (1998) [Prog. Part. Nucl. Phys. 41, 225 (1998)] 

arXiv:nucl-th/9803035. 
[14] M. Bleicher et ai, J. Phys. G 25, 1859 (1999) |arXiv:hep-ph/9909407 . 
[15] B. Nilsson-Almqvist and E. Stenlund, Comput. Phys. Commun. 43, 387 (1987). 
[16] T. Sjostrand, Comput. Phys. Commun. 82, 74 (1994). 
[17] M. Bleicher et ai, Nucl. Phys. A 638, 391 (1998). 

[18] F. Grassi, Y. Kama, O. Socolowski and T. Kodama, J. Phys. G 31, S1041 (2005). 

[19] R. Andrade, F. Grassi, Y. Kama, T. Kodama, O. . J. Socolowski and B. Tavares, Eur. Phys. J. A 

29, 23 (2006) [arXiv: nucl-th /051102l . 
[20] R. Andrade, F. Grassi, Y. Kama, T. Kodama and O. . J. Socolowski, Phys. Rev. Lett. 97, 202302 

(2006) arXiv:nucl-th/0608067 . 
[21] B. M. Tavares, H. J. Drescher and T. Kodama, Braz. J. Phys. 37, 41 (2007) 

(arXiv:hep-ph/0702224 . 

[22] R. P. G. Andrade, F. Grassi, Y. Kama, T. Kodama and W. L. Qian, Phys. Rev. Lett. 101, 112301 

(2008) arXiv:0805.0018 [hep-ph]]. 
[23] H. Petersen, T. Renk and S. A. Bass. larXiv: 1008.38461 [nucl-th]. 

[24] D. H. Rischke, S. Bernard and J. A. Maruhn, Nucl. Phys. A 595, 346 (1995) 
arXiv:nucl-th/95040l8 . 

[25] D. H. Rischke, Y. Pursun and J. A. Maruhn, Nucl. Phys. A 595, 383 (1995) [Erratum-ibid. A 596, 

717 (1996)] |arXiv:nucl-th/950402 1|. 
[26] H. Petersen, J. Steinheimer, M. Bleicher and H. Stocker, J. Phys. G 36, 055104 (2009) 

;arXiv:0902.4866' [nucl-th]] . 
[27] Anthony O'Hagan, Reliability Engineering k System Safety, 91 (10-11):1290-1300, (2006), The 

Fourth International Conference on Sensitivity Analysis of Model Output (SAMO 2004) - SAMO 

2004. 

[28] Jeremy E. Oakley and Anthony O'Hagan, Biometrika, 89 (4) :169-184, (2002). 
[29] Jeremy E. Oakley and Anthony O'Hagan, J. Roy. Stat. Soc. B,66,(2004). 

[30] M. J. Bayarri, James O. Berger, Rui Paulo, Jerome Sacks, John A. Cafeo, James C. Cavendish, 

Chin-Hsu Lin and Jian Tu, Technometrics, 49 (2):138-154, (2007). 
[31] David Higdon, James Gattiker, Brian Williams and Maria Rightley, J. Am. Stat. Assoc., 103 

(482):570-583,(2008). 
[32] Jean-Paul Chiles and Pierre Delfiner, John Wiley k Sons, New York, NY, (1999). 
[33] Noel A. C. Cressie, John Wiley & Sons, New York, NY, (1993). 
[34] Marc C. Kennedy and Anthony O'Hagan, Biometrika, 87(1):1-13, (2000). 

[35] Jerome Sacks, Wilham J. Welch, Toby J. Mitchell and Henry P. Wynn, Stat. Sci., 4 (4):409-435, 
November (1989). 

[36] Thomas J. Santner, Brian J. Williams and William Notz, Springer Series in Statistics, Springer- 

Verlag, New York, NY, (2003). 
[37] J. Adams et al. [STAR Collaboration], Phys. Rev. Lett. 92, 112301 (2004) arXivmucl-ex/0310004 . 
[38] S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. C 69, 034909 (2004) 
arXiv:nucl-ex/0307022; . 

[39] I. Arsene et al. [BRAHMS Collaboration], Phys. Rev. C 72, 014908 (2005) Jar'Xivmucl-ex/OSOMT^ . 
[40] S. Esumi [PHENIX CoUaboration] , Nucl. Phys. A 715, 599 (2003) arXivmucl-ex/021001 ^. 
[41] S. Manly et al. [PHOBOS Collaboration], Nucl. Phys. A 715, 611 (2003) arXiv:nucl-ex/ 0210036| . 
[42] R. L. Ray [STAR Collaboration], Nucl. Phys. A 715, 45 (2003) arXiv:nucl-ex/0211030 . 
[43] H. Petersen and M. Bleicher, Phys. Rev. C 81, 044906 (2010) iarXiv:1002.1003l [nucl-th]]. 



Initial State Granularity 



14 



[44] B. I. Abelev et at [STAR Collaboration], Phys. Rev. C 77, 054901 (2008) |arXiv:0801 .34661 [nucl- 

ex]]. ^ 

[45] H. Holopainen, H. Niemi and K. J. Eskola. larXiv: 1007.03681 [hep-ph]. 

[46] J. Steinheimer, V. Dexheimer, H. Petersen, M. Bleicher, S. Schramm and H. Stoecker, Phys. Rev. 

C 81, 044913 (2010) [arXiv:0905.3099 [hep-ph]]. 
[47] J. Steinheimer, S. Schramm and H. Stocker, arXiv:0909.4421 [hep-ph]. 

[48] S. Borsanyi, G. Endrodi, Z. Fodor et ai, JHEP 1011 (2010) 077 [arX iv:1007.2580l [hep-lat]]. 
[49] D. Ouerdane [BRAHMS CoUaboration] , Nucl. Phys. A 715, 478 (2003) |arXiv:nucl-ex/0212001] . 
[50] J. H. Lee et al. [BRAHMS Collaboration], J. Phys. G 30, S85 (2004). 

[51] J. Adams et al. [STAR CoUaboration], Phys. Rev. Lett. 98, 062301 (2007) |arXiv:nucl-ex/0606014| . 



